Vibrational properties of amorphous silicon from tight-binding 0(N) calculation 
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We present an 0{N) algorithm to study the vibrational properties of amorphous silicon within 
the framework of tight-binding approach. The dynamical matrix elements have been evaluated 
numerically in the harmonic approximation exploiting the short-range nature of the density matrix 
to calculate the vibrational density of states which is then compared with the same obtained from 
a standard 0(A'^*) algorithm. For the purpose of illustration, an 1000-atom model is studied to 
calculate the localization properties of the vibrational eigenstates using the participation numbers 
calculation. 
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In recent years there has been a tremendous progress in 
the development of efficient 0{N) methods for electronic 
structure calculation of materials. There exist a number 
of 0{N) methods which have been developed to calcu- 
late the total energy of the systems and the electronic 
forces ||, |[ ^, H, ^. Most of these algorithms utiHze, 
in one way or another, the fact that the density matrix of 
the system to be studied is short-ranged 0, 1). As a re- 
sult, calculation of local electronic properties are possible 
by introducing a localization volume surrounding the re- 
gion of interest where the local electronic properties are 
to be studied. The size of the localization volume de- 
pends on the decaying nature of the density matrix and 
for systems with a clean band gap, e.g., semi-conductors 
it has been observed that the density matrix decays ex- 
ponentially H, 1^. Hence it is possible to define a lo- 
calization volume having a linear dimension of the order 
of localization length of the density matrix. Although 
this docs not appear to be trivially true in metallic sys- 
tems [0, the locality in electronic structure calculation 
does indeed exist irrespective of the nature of the sys- 
tem and was first utilized by Friedel who introduced 
the term local density of states. Later Heine |^] extended 
this idea of local electronic structure calculation of solids 
in real space in the context of Recursion method 
which has been further elaborated by Kohn |ll] as the 
principle of nearsightedness of equilibrium systems. 

In this paper we present an efficient algorithm for the 
calculation of vibrational properties of amorphous silicon 
from a tight-binding p^ , [T^ Hamiltonian exploiting the 
short-range nature of the density matrix. The motivation 
behind this work is the calculation of localization prop- 
erties of the vibrational eigenstates of amorphous silicon. 
Recently, there has been a renewed interest in the study 
of the vibrational modes in disordered media jl^, ^ . 
These modes can be classified as either extended or local- 
ized depending on the energy of the state. In d = 3, the 
modes are expected to separate by sharp energy bound- 
ary, known as 'mobility edge' which marks the 'localized 
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to delocalized' transition. The corresponding analog in 
the electronic system is well-known after Anderson's sug- 
gestion jl^ that, unlike crystalline systems where elec- 
tronic states are extended over the system, sufficient dis- 
order can localize electronic states within the allowed en- 
ergy bands. For vibrational eigenstates, however, the low 
frequency, long wavelength modes (sound) always exist. 
Garber et al. studied an amorphous model of 4096 sil- 
icon atoms and demonstrated that at the mobility edge 
the vibrational eigenstates decay exponentially and lo- 
calization length diverges as a power law above the edge. 
Although we are not going to address these issues here, 
the importance of the construction of dynamical matrix 
from tight-binding calculations and to study the nature 
of the vibrational states for a large, well-relaxed sample 
of amorphous silicon has its own merits. 

In the following, we illustrate our work in two steps: 
first, we calculate the approximate electronic forces using 
an 0{N) algorithm starting from a tight-binding Hamil- 
tonian. This is then followed by a calculation of harmonic 
normal modes of vibrations by direct diagonalization of 
the dynamical matrix. The diagonalization of a matrix 
is an expensive 0{N^) operation, but since the time re- 
quired to construct the dynamical matrix elements using 
this 0{N) approach is much larger compared to a sin- 
gle diagonalization for the system sizes to be considered 
here, we make no attempt to obtain the density of states 
using any 0{N) algorithm. In addition, we also study 
the nature of the vibrational eigenstates for which it is 
necessary to compute the eigenvectors of the dynamical 
matrix for participation numbers calculation. For the 
density of states calculation only, however, one can use 
any one of the 0{N) algorithms available in the literature 
which include Maximum Entropy approach by Drabold 
and Sanke y Eg] , Recursion method of Haydock, Heine 
and Kelly |1C(| and Kernel Polynomial approximation by 
Silver and Roder etc. For our present purpose, 

a single diagonalization is easily affordable for the cal- 
culation of spectral quantities whereas the construction 
of a dynamical matrix using the exact electronic forces 
by direct diagonalization is extremely inefficient and time 
consuming which scales as 0{N'^) making it a formidable 
task to go beyond 500 atoms. This direct diagonalization 
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approach is not practical and feasible for very large sys- 
tem size (> 10000 atoms). However, one can still obtain 
spectral density distribution using any one of the meth- 
ods mentioned earlier and a few eigen vectors in appro- 
priately chosen intervals using Lanzcos algorithm in an 
0{N) way to study the nature of the states. 

Throughout this work, we have used the structural 
model generated by Barkema and Mousseau as start- 
ing configurations which have been relaxed using the 
tight-binding Hamiltonian parameterized by Kwon et 
al. [||. 

The general form of the tight-binding Hamiltonian 
H can be written as : H = J2i L^i^Llh L){i, L\ + 

J2z,LJ2i',L'^hL;i',L'\i:L){i',L'\ and the corresponding 

total energy is E ^ 2^21""'' {"fklHl^k) + Ur + Eat N. Ur 
stands for repulsive potential for the ion-ion interaction 
which can be represented either by a sum of pair poten- 
tials or as in the embedded atom model approach ||2^ 
as a functional of atomic density. N is the total num- 
ber of atoms and Eat is a constant energy /atom. The 
electronic eigen function \^k) oi H can be expanded as 
a linear combination : \'i'k) = J2i i of basis func- 

tions {|</';,i)}. For amorphous silicon, the minimal basis 
functions consist of one s-electron and three p-electrons 
per atom. 

The dynamical matrix element between two atoms at 
sites i and j can be written as 124] 
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where M^, Mj are the atomic masses at the sites i and j 
respectively. Fi'^ ^ F-^-'^ are the forces acting on the atom 
at the site j in the /? (= x, y, z) direction due to displace- 
ment Axi^Q. of the atom at the site i in the a direction 
before and after the displacement respectively. A straight 
forward implementation |l^ of the above procedure 
for obtaining dynamical matrix elements between all the 
atoms require at least (3A^-|-1) direct diagonahzations for 
a system of N number of atoms. The total force acting 
on an atom (say i) can be written as, 
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Using Hellman-Feynman theorem pq ], the electronic 
part can be written as : 



k j,L j',L' 



k j-L j',L' ' 



The calculation of eigen vectors (b's) is an 0{N^) oper- 
ation and for the entire system (containing N atoms) the 




FIG. 1: Partitioning of the entire space for the calculation 
of dynamical matrix elements between the atom at the cen- 
ter of the sphere and those inside the small sphere of radius 
Rd- The large sphere with radius Rl corresponds to the lo- 
calization volume (denoted by region I) for the atom at the 
center contributing to the forces on the atoms inside the small 
sphere. 



scaling goes as O(iV^). For each dynamical matrix ele- 
ment calculation one has to calculate the electronic force 
before and after the displacement of the atom and this 
has to be repeated for each direction /3(= x,y,z). For 
N > 1000, it is computationally overkill to construct the 
full dynamical matrix for the system. The calculation, 
however, can be performed in an 0{N) way by introduc- 
ing a localization volume surrounding the atom on which 
the force is to be calculated as described below. 

To calculate the electronic force using an 0{N) ap- 
proach, we partition the entire system into regions I and 
H with the the atom as center as shown in figure I. The 
large sphere with radius Rl corresponds to the localiza- 
tion volume for the central atom at i. The dynamical 
matrix elements between the atom i at the center and 
those inside the small sphere of radius Rd are to be con- 
structed. Typically, Rd can be chosen as second neighbor 
distance of the central atom or more. By varying the ra- 
dius Rd {Rd < Rl), one can monitor the magnitude of 
the dynamical matrix elements and can check whether 
they affect the vibrational density of states or not. The 
far away the atoms are from the central atom, the less is 
the effect of displacement on the atom at j [g6| . 

One can, in fact, introduce a cut-off for dynamical ma- 
trix elements while storing the elements. Here we have 
used a value of 3. 9 A for Rd and 10 A for Rl which we 
find is sufficient for the purpose of our calculation. The 
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FIG. 2: The convergence behavior of dynamical matrix ele- 
ments for different values of localization radius Rl- D and Do 
are the approximate and exact dynamical matrix elements re- 
spectively for each of the Rl values. For Rl > lOA, the error 
is observed to be less than 1 percent. 
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FIG. 3: The density of vibrational states for a system of 
500-atom model calculated using the 0(N) force algorithm 
described in the text. The same from a standard 0(A''^) force 
calculation is also plotted (denoted by a dotted line) for com- 
parison. 



difference of force on an atom j due to the displacement 
of the atom i at the center can be written as : 
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where F^oii^oii) ^^'^ li-^n ii) ^^'^ contribution 
to the total forces from the region I (region II) on the 
atom j before and after the displacement of the atom i 
at the center respectively. AFi and APi are the exact 
and approximate difference of total forces at the atom j. 
Our 0{N) algorithm is based on the assertion that the 
approximate difference of total force can be obtained by 
neglecting the contribution from region II if the radius 
of the region I is greater than the localization length of 
the density matrix. This assertion can be verified simply 
by looking at the convergence behavior of the dynam- 
ical matrix elements with respect to the radius of the 
localization volume. The approximate difference of force 
can now be written as a sum of electronic and repulsive 
contribution : 
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where Hf is the Hamiltonian of the region I before {x = 
o) and after the displacement (x = n) of the atom at the 
center. 

We have observed that a localization volume of radius 
lOA, typically containing 200 to 225 atoms, gives a very 
good estimate of APi for calculation of dynamical ma- 
trix elements. If we neglect contributions from region II, 
the force calculation on a single atom is an 0(1) oper- 
ation as the size of region I is independent of the total 
system size and the full dynamical matrix construction is 
an 0{N) operation for N atoms. Since the dimension of 
the truncated Hamiltonian matrices are much smaller 
than the Hamiltonian matrix H of the total system, the 
electronic forces can be calculated by direct diagonaliza- 
tion of the truncated Hamiltonian matrix Hj which is 
also independent of the system size. In practice, the ac- 
curacy of APi depends on two factors - the step size of 
displacement and the volume of region I which is again 
related with the size of band-gap. The accuracy can be 
improved in a controlled manner by (a) choosing a small 
displacement which is a necessary condition for the har- 
monic approximation to be valid and, (b) increasing the 
radius of the localization volume. We would like to men- 
tion that a similar but not identical approach based on 
the use of localized Wannier-like functions and truncat- 
ing the dynamical matrix beyond a cutoff was reported 
by Ordejon et al. [ p7[ . 

In figure |^ we have presented convergence behavior of 
dynamical matrix elements with respect to different lo- 
calization radius Rl. The estimate of approximate elec- 
tronic forces and hence the dynamical matrix elements 
are found to be very good and the percentage error is 
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FIG. 4: The vibrational density of states for a system of 1000- 
atom obtained using the 0{N) force algorithm described in 
the text. 



observed to be as low as 0.5 for Rl > lOA. The results 
of our calculations are presented in the figure ^ where we 
have plotted the density of vibrational states obtained 
from a 500-atom model. For the purpose of comparison, 
we have also plotted the density of states using the full 
Hamiltonian matrix by a standard 0{N'^) algorithm. It 
is quite clear from the figure || that the quality of (elec- 
tronic) forces obtained from the 0{N) algorithm is very 
good giving the same vibrational density of states as one 
would expect using the full Hamiltonian matrix. We 
should mention that to generate the dynamical matrix 
elements for this 500-atom model using standard 0{N^) 
force algorithm took 80 hours in contrast to 13.5 hours 
using the 0{N) algorithm in a DEC alpha (667 MHz) 
workstation. 

Having established that the dynamical matrix ob- 
tained in this 0{N) way is sufficiently accurate, we now 
study the nature of the vibrational states for a model of 
1000-atom. The configuration is first relaxed using the 
tight-binding Monte Carlo calculation as outlined in the 
Ref . . The vibrational density of states of the Monte 
Carlo relaxed structure is shown in the figure ^ which has 
been obtained by direct diagonalization of the dynami- 
cal matrix. Since the dynamical matrix construction for 
this 1000-atom model took approximately 27 hours, we 
therefore avoid to use any one of the 0{N) algorithm 
mentioned in the Refs. ||lj, |l9|, ^ for the calculation 
of vibrational density of states (VDOS). The VDOS in 
figure Q clearly shows two peaks : a low energy acousti- 
cal peak at 20 meV and a high energy optical peak at 80 
meV. The ratio of the height of two peaks are of the same 
order as observed in the experimental data obtained by 
Kamitakahara et al. The optical peak is, however, 

at higher energy compared to the experimental value and 
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FIG. 5: The participation numbers (in log scale) for a model 
of 1000-atom with energy. The number drops near around 80 
meV showing the localized character of the eigenstates which 
is marked by a vertical line. 



the value obtained from empirical potential calculations. 
Such a high energy optical peak has been also observed 
in the calculation of Klein et al. |^ who studied hydro- 
genated amorphous silicon using tight-binding molecular 
dynamics. We believe that this discrepancy in the po- 
sition of the high energy optical peak is associated with 
the repulsive part of the potential used in this calculation. 
The potential produces the correct electronic density of 
states within the accuracy of a fraction of an eV; how- 
ever, for the vibrational density of states calculation an 
accuracy of the order of few meVs is required to produce 
the peaks at right positions. 

In order to get an idea about the nature of the vi- 
brational states, we need to calculate the participation 
numbers from the eigenvectors of the dynamical matrix. 
Since the eigenvectors calculation is a memory as well as 
CPU intensive operation, it is not a very suitable tool 
for dealing with very large system sizes. We therefore re- 
strict ourselves to an 1000-atom model for this purpose. 
The participation numbers give a measure of the degree 
of localization of the wave functions to be studied. Let 
denotes the amplitude of the j-th normalized eigenvector. 
A simple measure of the number of sites participated in 
contributing the eigenstate gives the participation num- 
ber PnU) and is defined as P^^ = |0j(?)|^. A com- 
plete localized state therefore corresponds to Pn ~ 1, 
whereas for an extended state P/v = N. In figure |^, we 
have plotted the participation numbers with energy in 
meV. In the low energy region the value of the participa- 
tion number fluctuates around 300. The drop of the par- 
ticipation numbers around 80 meV is actually somewhat 
higher than that was observed in the work of Allen et 
al. who found a value of around 70 meV. This differ- 
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ence is, however, consistent with what we have observed 
in figure ^ where the optical peak appears at somewhat 
higher energy. In spite of this difference, we find that 
the general trend of the participation numbers with en- 
ergy (meV) follows the same course as was observed in 
Rcf. [|l4|. For example, a small dip at 30 mcV in the fig- 
ure ||has been also observed in Ref. |Q. A closer inspec- 
tion of the localization properties using the participation 
numbers requires a much larger system size to be studied 
for better statistics, particularly, in the high energy re- 
gion where a delocalized - localized transition is expected 
to occur from the low energy side. The tight-binding cal- 
culation for a model of 10 000-atom is in progress and in 
a future communication we will address the energy-level- 
spacing distribution and eigen function statistics of the 
spectrum. 

In conclusion, we have presented an 0{N) algorithm 
for the calculation electronic force for amorphous semi- 
conductors from tight-binding Hamiltonian exploiting 
the short-range nature of the density matrix. This algo- 
rithm has been used to construct the dynamical matrix 



for a 500-atom model of amorphous silicon. The vibra- 
tional density of states calculated from this 0{N) dynam- 
ical matrix has been observed to be in excellent agree- 
ment with the same obtained from a standard 0{N'^) 
algorithm. Furthermore, we have used this algorithm to 
study the nature of the vibrational states for a model 
of 1000-atom amorphous silicon, the dynamical matrix 
of which is otherwise difficult to construct in the tight- 
binding approach. The participation numbers calculation 
suggests that the results of this tight-binding calcula- 
tion are in agreement with the results based on empirical 
model potentials. 
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